.ipynb .pdf

repository open issue

Binder

Contents

Heat Transport

We will assume the heat budget in the model takes the form

(5)\[\begin{align} \frac{D \Theta}{D t} = - \nabla \cdot (\mathbf{F}^\Theta_{diff}) - \frac{1}{\rho_0 c_p }\frac{\partial Q}{\partial z} \end{align}\]

The term on the LHS represents the advective flux of conservative temperature. The first term on the RHS is the convergence of diffusive fluxes of heat. In the model, this represents turbulent diffusion due to mesoscale eddies and small-scale mixing. The final term is the convergence of heat fluxes from the atmosphere.

We will further assume a rigid lid. Note that the ECCO model actually uses a sophisticated nonlinear free surface form of the conservation equations, so these approximations aren’t quite right. But our goal here is to focus on the big picture, and these details would be a distraction.

import xgcm
import xarray as xr
from matplotlib import pyplot as plt
import numpy as np
import hvplot.xarray

from ecco_utils import load_ecco_v4

plt.rcParams['figure.figsize'] = (12, 6)

ds = load_ecco_v4()
ds

<xarray.Dataset>
Dimensions:    (face: 13, lon_c: 360, lon_g: 360, lat_c: 240, lat_g: 240, Z: 50, Zl: 50, Zp1: 51, Zu: 50, time: 288, time_snp: 287)
Coordinates: (12/42)
  * face       (face) int64 0 1 2 3 4 5 6 7 8 9 10 11 12
    i          (lon_c) int64 0 1 2 3 4 5 6 7 ... 352 353 354 355 356 357 358 359
    i_g        (lon_g) int64 0 1 2 3 4 5 6 7 ... 352 353 354 355 356 357 358 359
    j          (lat_c) int64 30 31 32 33 34 35 36 ... 264 265 266 267 268 269
    j_g        (lat_g) int64 30 31 32 33 34 35 36 ... 264 265 266 267 268 269
    k          (Z) int64 0 1 2 3 4 5 6 7 8 9 ... 40 41 42 43 44 45 46 47 48 49
    ...         ...
    dyG        (lat_c, lon_g) float32 dask.array<chunksize=(60, 90), meta=np.ndarray>
    dxG        (lat_g, lon_c) float32 dask.array<chunksize=(60, 90), meta=np.ndarray>
  * lon_c      (lon_c) float32 -37.5 -36.5 -35.5 -34.5 ... 319.5 320.5 321.5
  * lon_g      (lon_g) float32 -38.0 -37.0 -36.0 -35.0 ... 319.0 320.0 321.0
  * lat_c      (lat_c) float32 -81.46 -81.13 -80.8 -80.45 ... 67.16 67.34 67.47
  * lat_g      (lat_g) float32 -81.61 -81.28 -80.95 -80.61 ... 67.06 67.25 67.4
Data variables: (12/35)
    ADVr_SLT   (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    ADVr_TH    (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    DFrE_SLT   (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    DFrE_TH    (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    DFrI_SLT   (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    DFrI_TH    (time, Zl, lat_c, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    ...         ...
    UVELMASS   (time, Z, lat_c, lon_g) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    VVELMASS   (time, Z, lat_g, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    UVELSTAR   (time, Z, lat_c, lon_g) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    VVELSTAR   (time, Z, lat_g, lon_c) float32 dask.array<chunksize=(1, 50, 60, 90), meta=np.ndarray>
    oceTAUX    (time, lat_c, lon_g) float32 dask.array<chunksize=(1, 60, 90), meta=np.ndarray>
    oceTAUY    (time, lat_g, lon_c) float32 dask.array<chunksize=(1, 60, 90), meta=np.ndarray>

Vertically Integrated Heat Budget

Integrating the budget vertically and assuming \(\eta = 0\), we obtain

\[ \frac{\partial}{\partial t} \int_{-H}^0 \Theta dz = - \frac{\partial}{\partial x} \int_{-H}^0 u \Theta dz - \frac{\partial}{\partial y} \int_{-H}^0 v \Theta dz - \frac{\partial}{\partial x} \int_{-H}^0 F^\Theta_{diff,x} dz - \frac{\partial}{\partial y} \int_{-H}^0 F^\Theta_{diff,y} dz + \frac{Q}{\rho_0 c_p} \]

We can write this in a more compact form as $$ \frac{\partial}{\partial t} \int_{-H}^0 \Theta dz =

• \nabla_H \cdot (\mathbf{F}^\Theta_{adv,H} + \mathbf{F}^\Theta_{diff,H})

• \frac{Q}{\rho_0 c_p} $$

Where \(\mathbf{F}^\Theta_{adv,H}\) and \(\mathbf{F}^\Theta_{diff,H}\) are the vertically integrated advective and diffusive (2D) heat flux vectors. The operator \(\nabla_H\) represents the divergence in two dimensions (x, y).

from dask_gateway import Gateway
gway = Gateway()
cluster = gway.new_cluster()
cluster

grid = xgcm.Grid(ds, periodic=['X'])

fluxes = xr.Dataset({
    'F_adv_x': ds.ADVx_TH.sum('Z'),
    'F_adv_y': ds.ADVy_TH.sum('Z'),
    'F_diff_x': ds.DFxE_TH.sum('Z'),
    'F_diff_y': ds.DFyE_TH.sum('Z'),    
}).reset_coords(drop=True)
fluxes

<xarray.Dataset>
Dimensions:   (lon_g: 360, lat_c: 240, time: 288, lon_c: 360, lat_g: 240)
Coordinates:
  * time      (time) datetime64[ns] 1992-01-15 1992-02-13 ... 2015-12-14
  * lon_g     (lon_g) float32 -38.0 -37.0 -36.0 -35.0 ... 319.0 320.0 321.0
  * lat_c     (lat_c) float32 -81.46 -81.13 -80.8 -80.45 ... 67.16 67.34 67.47
  * lon_c     (lon_c) float32 -37.5 -36.5 -35.5 -34.5 ... 319.5 320.5 321.5
  * lat_g     (lat_g) float32 -81.61 -81.28 -80.95 -80.61 ... 67.06 67.25 67.4
Data variables:
    F_adv_x   (time, lat_c, lon_g) float32 dask.array<chunksize=(1, 60, 90), meta=np.ndarray>
    F_adv_y   (time, lat_g, lon_c) float32 dask.array<chunksize=(1, 60, 90), meta=np.ndarray>
    F_diff_x  (time, lat_c, lon_g) float32 dask.array<chunksize=(1, 60, 90), meta=np.ndarray>
    F_diff_y  (time, lat_g, lon_c) float32 dask.array<chunksize=(1, 60, 90), meta=np.ndarray>

try:
    fluxes = xr.open_dataset('tmp_data/ECCO_lateral_heat_fluxes.nc')
except FileNotFoundError:
    import dask
    cluster.scale(20)
    with cluster.get_client():
        for vname in fluxes.data_vars:
            fluxes[vname].load()
    cluster.scale(0)
    fluxes.to_netcdf('tmp_data/ECCO_lateral_heat_fluxes.nc')
fluxes

<xarray.Dataset>
Dimensions:   (time: 288, lon_g: 360, lat_c: 240, lon_c: 360, lat_g: 240)
Coordinates:
  * time      (time) datetime64[ns] 1992-01-15 1992-02-13 ... 2015-12-14
  * lon_g     (lon_g) float32 -38.0 -37.0 -36.0 -35.0 ... 319.0 320.0 321.0
  * lat_c     (lat_c) float32 -81.46 -81.13 -80.8 -80.45 ... 67.16 67.34 67.47
  * lon_c     (lon_c) float32 -37.5 -36.5 -35.5 -34.5 ... 319.5 320.5 321.5
  * lat_g     (lat_g) float32 -81.61 -81.28 -80.95 -80.61 ... 67.06 67.25 67.4
Data variables:
    F_adv_x   (time, lat_c, lon_g) float32 ...
    F_adv_y   (time, lat_g, lon_c) float32 ...
    F_diff_x  (time, lat_c, lon_g) float32 ...
    F_diff_y  (time, lat_g, lon_c) float32 ...

over_rA = (ds.rA.compute()**-1).fillna(0.)

advective_flux_convergence = -(
    grid.diff(fluxes.F_adv_x, 'X', boundary='extend') +
    grid.diff(fluxes.F_adv_y, 'Y', boundary='extend')
) * over_rA

diffusive_flux_convergence = -(
    grid.diff(fluxes.F_diff_x, 'X', boundary='extend') +
    grid.diff(fluxes.F_diff_y, 'Y', boundary='extend')
) * over_rA

rho0 = 1029.
cp = 3850.

land_mask = (ds.hFacC[0] > 0).compute()
surface_flux = ds.TFLUX.compute().where(land_mask) / (rho0 * cp)

clim = (-1e-4, 1e-4)
big_plot = (
 advective_flux_convergence.hvplot.quadmesh('lon_c', 'lat_c', clim=clim, cmap='RdBu_r', label='Advective_Flux_Convergence') +
 diffusive_flux_convergence.hvplot.quadmesh('lon_c', 'lat_c', clim=clim, cmap='RdBu_r', label='Meridional_Flux_Convergence') +
 surface_flux.hvplot.quadmesh('lon_c', 'lat_c',clim=clim, cmap='RdBu_r', label='Surface_Flux') 
).cols(1)
big_plot

In the long term mean, we see that the convergence of the lateral fluxes balances the net air-sea fluxes quite closely.

import cartopy.crs as ccrs
proj = ccrs.Robinson(central_longitude=160)

fig, axes = plt.subplots(ncols=2, figsize=(16, 5), subplot_kw=dict(projection=proj, facecolor='0.9'))

for ax in axes:                                                       
    ax.coastlines()

clim = (-2e-5, 2e-5)

pc1 = surface_flux.mean('time').plot(ax=axes[0], transform=ccrs.PlateCarree(), add_colorbar=False)
pc2 = (advective_flux_convergence + diffusive_flux_convergence).mean('time').plot(ax=axes[1], transform=ccrs.PlateCarree(), add_colorbar=False)
pc1.set_clim(clim)
pc2.set_clim(clim)
axes[0].set_title('Time Mean Surface Flux')
axes[1].set_title('Time Mean Heat Transport Convergence')
plt.close()

fig

Ocean Meridional Heat Transport

If we integrate one more time in the \(x\) (i.e. Longitude) direction, around closed latitude circules we obtain a heat budget that is a function of latitude and time only.

\[ \frac{\partial}{\partial t} \oint dx \int_{-H}^0 \Theta dz = - \frac{\partial}{\partial y} \oint \int_{-H}^0 ( v \Theta + \mathbf{F}^\Theta_{diff, x}) dz dx + \oint \frac{Q}{\rho_0 c_p} dx \]

The first term RHS is the Ocean Meridional Heat Transport, a fundamental quantity in the climate system.

F_merid_adv = (rho0 * cp * fluxes.F_adv_y.sum('lon_c') / 1e15)
F_merid_diff = (rho0 * cp * fluxes.F_diff_y.sum('lon_c') / 1e15)
F_merid_tot = (F_merid_adv + F_merid_diff)

F_merid_tot.mean('time', keep_attrs=True).plot(color='black', label='Meridional Heat Transport')
F_merid_adv.mean('time', keep_attrs=True).plot(label='Advective Meridional Heat Transport')
F_merid_diff.mean('time', keep_attrs=True).plot(label='Diffusive Meridional Heat Transport')
plt.legend()
plt.grid()
plt.ylabel('PW')

Text(0, 0.5, 'PW')

Separating Atlatic and Indo-Pacific

atl_basins = [2, 7, 9, 10, 12, 19]
atl_mask = sum([ds.basins == basin_num for basin_num in atl_basins])

mask_c = atl_mask.compute()
mask_v = (grid.interp(mask_c, 'Y', boundary='fill') > 0)

F_merid_atl = rho0 * cp * (fluxes.F_diff_y + fluxes.F_adv_y).where(mask_v).sum('lon_c').sel(lat_g=slice(-30, None)) / 1e15
F_merid_rest = rho0 * cp * (fluxes.F_diff_y + fluxes.F_adv_y).where(~mask_v).sum('lon_c').sel(lat_g=slice(-30, None)) / 1e15

F_merid_tot.mean('time').plot(color='black', label='Meridional Heat Transport')
F_merid_adv.mean('time').plot(label='Advective Meridional Heat Transport')
F_merid_atl.mean('time').plot(label='Atlantic')
F_merid_rest.mean('time').plot(label='Indo Pacific')
plt.legend()
plt.grid()
plt.ylabel('PW')

Text(0, 0.5, 'PW')

Decomposing the Advective Meridional Heat Transport

(6)\[\begin{align} \mathcal{F}_{adv} &= \oint \int_{-H}^0 v \Theta dz dx \\ \end{align}\]

Now let’s define the zonal mean velocity and temperature as

\[ \overline{v} = \frac{1}{L} \oint v dx \ , \ \ \overline{\Theta} = \frac{1}{L} \oint \Theta dx\]

where \(L\) is the length of the latitude circle.

The total velocity and temperature can be written as a sum of this mean component and a fluctuation from it.

\[ v = \overline{v} + v' \ , \ \ \Theta = \overline{\Theta} + \Theta' \]

So

(7)\[\begin{align} \mathcal{F}_{adv} &= L_x \int_{-H}^0 \overline{v} \overline{\Theta} dz dx + \oint \int_{-H}^0 \overline{v' \Theta'} dz dx \\ \end{align}\]

ds_transport = xr.open_dataset('tmp_data/ECCO_overturning_transport.nc')
ds_transport

<xarray.Dataset>
Dimensions:      (time: 288, Z: 50, lat_g: 240)
Coordinates:
  * time         (time) datetime64[ns] 1992-01-15 1992-02-13 ... 2015-12-14
  * Z            (Z) float32 -5.0 -15.0 -25.0 ... -5.461e+03 -5.906e+03
  * lat_g        (lat_g) float32 -81.61 -81.28 -80.95 ... 67.06 67.25 67.4
Data variables:
    V_atl        (time, Z, lat_g) float32 ...
    V_rest       (time, Z, lat_g) float32 ...
    V_eddy_atl   (time, Z, lat_g) float32 ...
    V_eddy_rest  (time, Z, lat_g) float32 ...

with cluster.get_client():
    Theta_bar = ds.THETA.where(ds.hFacC>0).mean('lon_c').reset_coords(drop=True).compute()
    Theta_bar_atl = ds.THETA.where(mask_c).where(ds.hFacC>0).mean('lon_c').reset_coords(drop=True).compute()
    Theta_bar_rest = ds.THETA.where(~mask_c).where(ds.hFacC>0).mean('lon_c').reset_coords(drop=True).compute()

Theta_bar_g = grid.interp(Theta_bar, 'Y', boundary='extend')
Theta_bar_atl_g = grid.interp(Theta_bar_atl, 'Y', boundary='extend')
Theta_bar_rest_g = grid.interp(Theta_bar_rest, 'Y', boundary='extend')

F_adv_MOC = rho0 * cp * (Theta_bar_g * (ds_transport.V_atl + ds_transport.V_rest)).sum('Z') / 1e15
F_adv_MOC_atl = rho0 * cp * (Theta_bar_atl_g * ds_transport.V_atl).sum('Z').sel(lat_g=slice(-30, None)) / 1e15
F_adv_MOC_rest = rho0 * cp * (Theta_bar_rest_g * ds_transport.V_rest).sum('Z').sel(lat_g=slice(-30, None)) / 1e15

F_adv_MOC.mean('time').plot()
F_adv_MOC_atl.mean('time').plot()
F_adv_MOC_rest.mean('time').plot()
(F_adv_MOC_atl + F_adv_MOC_rest).mean('time').plot()

[<matplotlib.lines.Line2D at 0x7ff6922efb80>]

F_merid_tot.mean('time').plot(color='black', label='Meridional Heat Transport')
F_merid_adv.mean('time').plot(label='Advective Meridional Heat Transport', color='grey')
F_merid_atl.mean('time').plot(label='Atlantic', color='r')
F_merid_rest.mean('time').plot(label='Indo Pacific', color='b')
F_adv_MOC_atl.mean('time').plot(label='Atlantic MOC', color='r', linestyle=':')
F_adv_MOC_rest.mean('time').plot(label='Indo Pacific MOC', color='b', linestyle=':')
plt.legend()
plt.grid()
plt.ylabel('PW')

Text(0, 0.5, 'PW')

previous

Overturning Circulation